Effect of basis degree, penalty and knots

devtools::install("../../elasticproc2d", quiet=TRUE, quick=TRUE)

library(elasticproc2d)
library(ggplot2)
library(grid)
library(gridExtra)
library(dplyr)
library(cowplot)

# Display plots in smaller size
options(repr.plot.width=15, repr.plot.height=5)

Load some datasets with random rotation and scaling.

source("../scripts/load_datasets.R")

set.seed(18)

spirals <- curves.spiral(n_curves=10, rotate=TRUE, scale=TRUE, center=TRUE)
digits3 <- curves.digit3(rotate=TRUE, scale=TRUE, center=TRUE)

Normalize curves.

Take a look at the data curves.

p1 <- ggplot(bind_rows(spirals, .id="id"), aes(x=X1, y=X2)) + 
    geom_path(size=0.5, aes(group=id), color="grey") +
    coord_fixed()

p2 <- ggplot(bind_rows(digits3, .id="id"), aes(x=X1, y=X2)) + 
    geom_path(size=0.5, aes(group=id), color="grey") +
    coord_fixed()

grid.arrange(p1, p2, nrow=1, top="Spirals and digits3 with random rotation and scaling")

Define functions for plotting.

plot_cov_surface<- function(cov_fit, h=0.01){
  # Define covariance surface grid (s,t).
  arg.grid = seq(0, 1, by=h)
  cov.grid = expand.grid(t = arg.grid, s = arg.grid)
  # Evaluate fit on grid.
  cov.re = predict(cov_fit$re, newdata = cov.grid)
  cov.im = predict(cov_fit$im, newdata = cov.grid)

  par(mfrow=c(1,2), mar=c(4,4,4,1), oma=c(0.5,0.5,0.5,0))
  # From 'fdapace/src/R/CreateCovPlot.R'
  args1 <- list(
    xlab='t', ylab='s', zlab = 'C(s,t)',
    main = 'Smoothed covariance surface (real part)'
  )
  args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.re, nrow=101))
  do.call(plot3D::persp3D, c(args2, args1))

  # From 'fdapace/src/R/CreateCovPlot.R'
  args1 <- list(
    xlab='t', ylab='s', zlab = 'C(s,t)',

    main = 'Smoothed covariance surface (imaginary part)'
  )
  args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.im,nrow=101))
  do.call(plot3D::persp3D, c(args2, args1))
}

library(grid)


plot_all <- function(mean, h=0.01){
    
  cov_fit <- mean$fit$cov_fit
  # Define covariance surface grid (s,t).
  arg.grid = seq(0, 1, by=h)
  cov.grid = expand.grid(t = arg.grid, s = arg.grid)
  # Evaluate fit on grid.
  cov.re = predict(cov_fit$re, newdata = cov.grid)
  cov.im = predict(cov_fit$im, newdata = cov.grid)

  par(mfrow=c(1,3), mar=c(4,4,4,1), oma=c(0.33,0.33,0.33,0))
 
  # From 'fdapace/src/R/CreateCovPlot.R'
  args1 <- list(
    xlab='t', ylab='s', zlab = 'C(s,t)',
    main = 'Smoothed covariance surface (real part)'
  )
  args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.re, nrow=101))
  do.call(plot3D::persp3D, c(args2, args1))

  # From 'fdapace/src/R/CreateCovPlot.R'
  args1 <- list(
    xlab='t', ylab='s', zlab = 'C(s,t)',

    main = 'Smoothed covariance surface (imaginary part)'
  )
  args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.im,nrow=101))
  do.call(plot3D::persp3D, c(args2, args1))

  plot(mean)
}

Testing Parameters

order <- list("polygon", "smooth")
penalty <- list(0, 1, 2)

# Turn of warnings
defaultW <- getOption("warn")
options(warn = -1)

Spirals

spiral.means <- lapply(order, function(type){
    lapply(penalty, function(pen){
        start_time <- Sys.time()
        knots <- seq(0, 1, length = if_else(type == "smooth", 13, 19))
        print(paste0("Computing spirals mean:  Type = ", type, "  Penalty = ", pen, "  #Knots = ", length(knots)))
        flush.console()
        mean <- compute_elastic_proc2d_mean(spirals, knots = knots, type = type, penalty = pen)
        print(Sys.time() - start_time)
        flush.console()
        mean
    })
})

[1] "Computing spirals mean:  Type = polygon  Penalty = 0  #Knots = 19"
Time difference of 1.163191 mins
[1] "Computing spirals mean:  Type = polygon  Penalty = 1  #Knots = 19"
Time difference of 44.74742 secs
[1] "Computing spirals mean:  Type = polygon  Penalty = 2  #Knots = 19"
Time difference of 20.94888 secs
[1] "Computing spirals mean:  Type = smooth  Penalty = 0  #Knots = 13"
Time difference of 11.57825 secs
[1] "Computing spirals mean:  Type = smooth  Penalty = 1  #Knots = 13"
Time difference of 11.86377 secs
[1] "Computing spirals mean:  Type = smooth  Penalty = 2  #Knots = 13"
Time difference of 11.56085 secs

void <- lapply(1:length(order), function(i){
    lapply(1:length(penalty), function(j){
        mean <- spiral.means[[i]][[j]]
        plot_all(mean)
        print(paste0("Type = ", mean$type, "  Penalty = ", penalty[[j]], "  #Knots = ", length(mean$knots)))
    })
})

[1] "Type = polygon  Penalty = 0  #Knots = 19"

[1] "Type = polygon  Penalty = 1  #Knots = 19"

[1] "Type = polygon  Penalty = 2  #Knots = 19"

[1] "Type = smooth  Penalty = 0  #Knots = 13"

[1] "Type = smooth  Penalty = 1  #Knots = 13"

[1] "Type = smooth  Penalty = 2  #Knots = 13"

Digits3

digits3.means <- lapply(order, function(type){
    lapply(penalty, function(pen){
        start_time <- Sys.time()
        knots <- seq(0, 1, length = if_else(type == "smooth", 13, 19))
        print(paste0("Computing digits3 mean:  Type = ", type, "  Penalty = ", pen, "  #Knots = ", length(knots)))
        flush.console()
        mean <- compute_elastic_proc2d_mean(digits3, knots = knots, type = type, penalty = pen)
        print(Sys.time() - start_time)
        flush.console()
        mean
    })
})

[1] "Computing digits3 mean:  Type = polygon  Penalty = 0  #Knots = 19"
Time difference of 56.27127 secs
[1] "Computing digits3 mean:  Type = polygon  Penalty = 1  #Knots = 19"
Time difference of 20.34927 secs
[1] "Computing digits3 mean:  Type = polygon  Penalty = 2  #Knots = 19"
Time difference of 20.33894 secs
[1] "Computing digits3 mean:  Type = smooth  Penalty = 0  #Knots = 13"
Time difference of 32.41793 secs
[1] "Computing digits3 mean:  Type = smooth  Penalty = 1  #Knots = 13"
Time difference of 32.61988 secs
[1] "Computing digits3 mean:  Type = smooth  Penalty = 2  #Knots = 13"
Time difference of 32.37099 secs

void <- lapply(1:length(order), function(i){
    lapply(1:length(penalty), function(j){
        mean <- digits3.means[[i]][[j]]
        plot_all(mean)
        print(paste0("Type = ", mean$type, "  Penalty = ", penalty[[j]], "  #Knots = ", length(mean$knots)))
    })
})

[1] "Type = polygon  Penalty = 0  #Knots = 19"

[1] "Type = polygon  Penalty = 1  #Knots = 19"

[1] "Type = polygon  Penalty = 2  #Knots = 19"

[1] "Type = smooth  Penalty = 0  #Knots = 13"

[1] "Type = smooth  Penalty = 1  #Knots = 13"

[1] "Type = smooth  Penalty = 2  #Knots = 13"

No penalty

mean.sp.nofx.poly <- compute_elastic_proc2d_mean(spirals, knots = seq(0,1,length=11), type = "polygon", penalty = -1)
mean.sp.nofx.smooth <- compute_elastic_proc2d_mean(spirals, knots = seq(0,1,length=11), type = "smooth", penalty = -1)

mean.d3.nofx.poly <- compute_elastic_proc2d_mean(digits3, knots = seq(0,1,length=11), type = "polygon", penalty = -1)
mean.d3.nofx.smooth <- compute_elastic_proc2d_mean(digits3, knots = seq(0,1,length=11), type = "smooth", penalty = -1)

Plots

plot_all(mean.sp.nofx.poly)
plot_all(mean.sp.nofx.smooth)

plot_all(mean.d3.nofx.poly)
plot_all(mean.d3.nofx.smooth)